#pragma once

#include <stdexcept>
#include <cmath>
#include <iostream>
#include <sstream>

namespace BR
{
	template<class Tp>
	struct Dual
	{
		typedef Tp value_type;
		value_type real;
		value_type imag;

		// constructor
		Dual(value_type const & r = value_type(), value_type const & i = value_type()) : real(r), imag(i) { }

		Dual(Dual const & rhs) : real(rhs.real), imag(rhs.imag) { }

		template<class Up> Dual(Dual<Up> const & rhs) : real(rhs.real), imag(rhs.imag) { }

		// value
		Dual & real(value_type const & val) { real = val; return *this; }
		Dual & imag(value_type const & val) { imag = val; return *this; }

		value_type real(void) const { return real; }
		value_type imag(void) const { return imag; }

		// assign
		Dual & assign(value_type const & first, value_type const & second) { real = first; imag = second; return *this; }
		Dual & assign(value_type const & rhs) { real = rhs; imag = Tp(); return *this; }

		template<class Up> Dual & assign(Dual<Up> const & rhs) { real = rhs.real; imag = rhs.imag; return *this; }

		Dual & operator=(value_type const & rhs) { return assign(rhs); }

		template<class Up> Dual & operator=(Dual<Up> const & rhs) { return assign(rhs); }

		// operation
		template<class Up> Dual & self_add(Dual<Up> const & rhs) { real += rhs.real; imag += rhs.imag; return *this; }
		template<class Up> Dual & self_sub(Dual<Up> const & rhs) { real -= rhs.real; imag -= rhs.imag; return *this; }
		template<class Up> Dual & self_mul(Dual<Up> const & rhs) { imag = real*rhs.imag + imag*rhs.real; real *= rhs.real; return *this; }
		template<class Up> Dual & self_div(Dual<Up> const & rhs) { imag = (imag*rhs.real - real*rhs.imag) / rhs.real; real /= rhs.real; return *this; }

		Dual & self_add(value_type const & rhs) { real += rhs; return *this; }
		Dual & self_sub(value_type const & rhs) { real -= rhs; return *this; }
		Dual & self_mul(value_type const & rhs) { real *= rhs; imag *= rhs; return *this; }
		Dual & self_div(value_type const & rhs) { real /= rhs; imag /= rhs; return *this; }

		template<class Up> Dual add(Dual<Up> const & rhs) const { Dual tmp(*this); return tmp.self_add(rhs); }
		template<class Up> Dual sub(Dual<Up> const & rhs) const { Dual tmp(*this); return tmp.self_sub(rhs); }
		template<class Up> Dual mul(Dual<Up> const & rhs) const { Dual tmp(*this); return tmp.self_mul(rhs); }
		template<class Up> Dual div(Dual<Up> const & rhs) const { Dual tmp(*this); return tmp.self_div(rhs); }

		Dual add(value_type const & rhs) const { Dual tmp(*this); return tmp.self_add(rhs); }
		Dual sub(value_type const & rhs) const { Dual tmp(*this); return tmp.self_sub(rhs); }
		Dual mul(value_type const & rhs) const { Dual tmp(*this); return tmp.self_mul(rhs); }
		Dual div(value_type const & rhs) const { Dual tmp(*this); return tmp.self_div(rhs); }

		Dual pos(void) const { return Dual( x.real,  x.imag); }
		Dual neg(void) const { return Dual(-x.real, -x.imag); }
		Dual con(void) const { return Dual( x.real, -x.imag); }

		bool eql(Dual const & rhs) const { return real==rhs.real && imag==rhs.imag; }
		bool neq(Dual const & rhs) const { return real!=rhs.real || imag!=rhs.imag; }

		bool eql(value_type const & rhs) const { return real==rhs && imag==Tp(); }
		bool neq(value_type const & rhs) const { return real!=rhs || imag!=Tp(); }

		bool eql(value_type const & first, value_type const & second) const { return real==first && imag==second; }
		bool neq(value_type const & first, value_type const & second) const { return real!=first || imag!=second; }

		// operator
		template<class Up> Dual & operator+=(Dual<Up> const & rhs) { return self_add(rhs); }
		template<class Up> Dual & operator-=(Dual<Up> const & rhs) { return self_sub(rhs); }
		template<class Up> Dual & operator*=(Dual<Up> const & rhs) { return self_mul(rhs); }
		template<class Up> Dual & operator/=(Dual<Up> const & rhs) { return self_div(rhs); }

		Dual & operator+=(value_type const & rhs) { return self_add(rhs); }
		Dual & operator-=(value_type const & rhs) { return self_sub(rhs); }
		Dual & operator*=(value_type const & rhs) { return self_mul(rhs); }
		Dual & operator/=(value_type const & rhs) { return self_div(rhs); }

		template<class Up> Dual operator+(Dual<Up> const & rhs) const { return add(rhs); }
		template<class Up> Dual operator-(Dual<Up> const & rhs) const { return sub(rhs); }
		template<class Up> Dual operator*(Dual<Up> const & rhs) const { return mul(rhs); }
		template<class Up> Dual operator/(Dual<Up> const & rhs) const { return div(rhs); }

		Dual operator+(void) const { return pos(); }
		Dual operator-(void) const { return neg(); }
		Dual operator~(void) const { return con(); }

		bool operator==(Dual const & rhs) const { return eql(rhs); }
		bool operator!=(Dual const & rhs) const { return neq(rhs); }

		value_type const & operator[](int idx) const
		{
			switch (idx) {
				case 0: return real;
				case 1: return imag;
				default: throw std::out_of_range("Invalid index.\n");
			}
		}

		value_type & operator[](int idx)
		{
			switch (idx) {
				case 0: return real;
				case 1: return imag;
				default: throw std::out_of_range("Invalid index.\n");
			}
		}
	};

	template<class Tp> inline Dual<Tp> operator+(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.add(rhs); }
	template<class Tp> inline Dual<Tp> operator-(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.sub(rhs); }
	template<class Tp> inline Dual<Tp> operator*(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.mul(rhs); }
	template<class Tp> inline Dual<Tp> operator/(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.div(rhs); }

	template<class Tp> inline Dual<Tp> operator+(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.add(lhs); }
	template<class Tp> inline Dual<Tp> operator-(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.sub(lhs); }
	template<class Tp> inline Dual<Tp> operator*(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.mul(lhs); }
	template<class Tp> inline Dual<Tp> operator/(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.div(lhs); }

	template<class Tp> inline bool operator==(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.eql(rhs); }
	template<class Tp> inline bool operator!=(Dual<Tp> const & lhs, Tp const & rhs) const { return lhs.neq(rhs); }

	template<class Tp> inline bool operator==(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.eql(lhs); }
	template<class Tp> inline bool operator!=(Tp const & lhs, Dual<Tp> const & rhs) const { return rhs.neq(lhs); }

	template<class Tp> inline Tp & real(Dual<Tp> & z) { return z.real; }
	template<class Tp> inline Tp & imag(Dual<Tp> & z) { return z.imag; }

	template<class Tp> inline Tp abs(Dual<Tp> const & z) { return abs(z.real); }

	template<class Tp> inline Tp arg(Dual<Tp> const & z) { return z.imag/z.real; }

	template<class Tp> inline Tp norm(Dual<Tp> const & z) { return z.real*z.real; }

	template<class Tp> inline Dual<Tp> conj(Dual<Tp> const & z) { return z.con; }

	template<class Tp> inline Dual<Tp> polar(Tp const & rho, Tp const & theta) { return Dual<Tp>(rho, rho*theta); }

	using std::sin;
	using std::cos;
	using std::tan;
	using std::sinh;
	using std::cosh;
	using std::tanh;
	using std::exp;
	using std::log;
	using std::pow;

	template<class Tp> inline Dual<Tp> sin (Dual<Tp> const & z) { return Dual<Tp>(sin (z.real),  z.imag * cos (z.real)); }
	template<class Tp> inline Dual<Tp> cos (Dual<Tp> const & z) { return Dual<Tp>(cos (z.real), -z.imag * sin (z.real)); }
	template<class Tp> inline Dual<Tp> sinh(Dual<Tp> const & z) { return Dual<Tp>(sinh(z.real),  z.imag * cosh(z.real)); }
	template<class Tp> inline Dual<Tp> cosh(Dual<Tp> const & z) { return Dual<Tp>(cosh(z.real),  z.imag * sinh(z.real)); }
	template<class Tp> inline Dual<Tp> exp (Dual<Tp> const & z) { return Dual<Tp>(exp (z.real),  z.imag * exp (z.real)); }
	template<class Tp> inline Dual<Tp> log (Dual<Tp> const & z) { return Dual<Tp>(log (z.real),  z.imag / z.real); }

	template<class Tp> inline Dual<Tp> tan (Dual<Tp> const & z) { Tp t(cos (z.real)); return Dual<Tp>(tan (z.real), z.imag/(t*t)); }
	template<class Tp> inline Dual<Tp> tanh(Dual<Tp> const & z) { Tp t(cosh(z.real)); return Dual<Tp>(tanh(z.real), z.imag/(t*t)); }

	template<class Tp> inline Dual<Tp> pow(Dual<Tp> const & z, int k) { return Dual<Tp>(pow(z.real, k), k*z.imag*pow(z.real, k-1)); }

	template<class Tp> inline Dual<Tp> log10(Dual<Tp> const & z)
	{
		static Tp const CONSTANT(log(Tp(10)));
		return Dual<Tp>(log (z.real),  z.imag / z.real / CONSTANT);
	}

	template<class Tp> inline Dual<Tp> sqrt(Dual<Tp> const & z)
	{
		static Tp const CONSTANT(2);
		return Dual<Tp>(sqrt(z.real),  z.imag / sqrt(z.real) / CONSTANT);
	}

	template<class Tp, class CharT, class Traits>
	std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits> &, Dual<Tp> &);

	template<class Tp, class CharT, class Traits>
	std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits> &, Dual<Tp> const &);
}
